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Abstract 

Approximate p-point Leibniz derivation formulas as well as interpolatory 
Simpson quadrature sums adapted to oscillatory functions are discussed. Both 
theoretical considerations and numerical evidence concerning the dependence of 
the discretization errors on the frequency parameter of the oscillatory functions 
show that the accuracy gain of the present formulas over those based on the 
exponential fitting approach [L. Ixaru, Computer Physics Communications, 105 
(1997) 1-19] is overwhelming. 



1 Introduction 

The mathematical descriptions of classical oscillatory phenomena (like vibrations, wave 
propagation, resonances) or of the behaviour of quantum systems involve operations 
on oscillatory functions. In most instances, such operations (e.g., differentiation, inte- 
gration, solving differential equations) require numerical methods. 

A successful approach towards accurate approximation of the oscillatory solutions 
of ordinary differential equations is the exponential fitting method, first proposed for 
the radial Schrodinger equation and then extended by several authors to more 
general differential equations (see, e.g., for a list of relevant results). 

* E-mail: adamg(8)theorl. theory.nipne.ro 
^ E-mail: adams@theorl.theory.nipne.ro 
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In a recent paper, Ixaru raised the question whether the exponential fitting 
approach could yield useful formulas for the numerical derivation and integration of 
oscillatory functions as well. His investigation resulted in several new formulas which 
extended to oscillatory functions well known elementary discrete approximations of the 
first and second order derivatives as well as the Simpson quadrature formula. 

With the aim at enhancing code robustness and reliability over integration sub- 
ranges characterized by the occurrence of a slowly varying regular factor, we tried to 
incorporate the optimized Simpson quadrature formula of in our computer code 
devoted to the integration of the products of functions with oscillatory factors 0, H, ^ 
within a class conscious automatic adaptive quadrature frame 0. Since the attempt 
resulted in significant deterioration of code performances, we undertook a complemen- 
tary study the results of which are reported below. 

There are four fundamental functional dependences of interest, namely, 

n^) = Kn^^) = fi^)9sA^^ + 5) , (1) 

where f{x) is a sufficiently smooth real function, while g^ ,^ {s = 1,2 , rj = ±1) denotes 
one of the following four weight functions 

gi^^i{ujx + 6) = cos{u!x + 6) , g2~i{ujx + 6) = sm{u!x + 6) , , . 

gi^i{ujx + 5) = cosh{ujx + 6) , g2,i{ujx + 6) = sinh(c<jx + 5) . 

Here, the frequency parameter u and the initial phase 6 are both real and constant. 

The principle of the present approach towards numerical differentiation and nu- 
merical integration of functions of the form ([^) or of their linear combinations is well- 
known (see, e.g., the derivation of extended Clenshaw-Curtis quadrature sums for oscil- 
latory functions in |^): the approximating operations apply to the regular factor f{x) 
only, while the contribution of the weight factor gs^rj is included exactly. This results in 
formulas which are uniformly valid at values of the frequency parameter u running over 
the whole set of machine numbers at which the pair of weight functions {5^1^^ (wx + 6), 
g2^r]{^x + S)} can be accurately computed. Of course, the classical differentiation and 
integration formulas are recovered in the limit uj —>■ 0. 

The derivation of exponential fitting formulas of numerical differentiation and 
integration of |0, requires the preservation of the formal structure of the corresponding 
classical formulas. The consequence is that the exponentially fitted coefficients depend 
on the gs^rf factor in a manner which results in breakdowns of the obtained formulas 
at some specific frequency parameter values. Due to this feature, the formulas of 
numerical differentiation and integration based on exponential fitting are of limited 
practical value. In particular, this is the reason for the abovementioned performance 
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deterioration of the code of automatic quadrature under incorporation of the optimal 
exponential fitting generalization of Simpson quadrature formula. 

To allow straightforward comparison of the uj dependences of the discretization 
errors associated to the approximating formulas within the two approaches, the case 
study problems considered in have been solved by both methods. The accuracy 
improvement brought by the present formulas with respect to those based on exponen- 
tial fitting was found to be comparable to that reported in to be brought by the 
exponential fitting ones with respect to their classical counterparts. 

The paper is organized in three sections. Section |^ discusses discretized Leibniz 
derivation formulas. Section ^ provides three-point interpolatory quadrature sums for 
integrands of the form (|l]). Section ^ summarizes the main practical consequences of 
the present investigation. 



2 Derivatives 



2.1 Discretized Leibniz formulas 

The description of oscillatory phenomena generally involves series of harmonics, with 
terms of the form 

<l>(x) = fi{x)gi^r^{ujx + 5) + f2{x)g2,n{'^x + 5) . (3) 

To get numerical differentiation or integration formulas of such an expression, solutions 
for the basic functional dependences (P are needed. 

Straightforward use of the Leibniz derivation formula yields for the n-th derivative 

of® 

(x) = E ( ^ ) ^'f^""'^ (^)^S (^) ' (4) 

where, for each involved function y{u)^ y^"^\u) = d^y{u) /du"^, while t = cox + 5 . For 
any natural number k, 

g^k+s)^t) = r^-i gs-sAt), 9i%'^'\t) = g,,,{t), 

hence (H) consists of a superposition of the pair of functions gi,r,{t) and g2,ri(t), with co- 
efficients which are themselves expressed as superpositions of derivatives of the regular 
factor f{x). 



3 



The use of the Leibniz formula (Q), combined with classical p-point difference 
classical discretization formulas for the involved derivatives of the regular factors /i (x) 
and f2ix) respectively will result in approximate derivation formulas of Eq. (H), 

which are optimal with respect to the frequency parameter u. The p-point discretiza- 
tions of the first and second order derivatives of are considered in detail below. 



2.1.1 First Order Derivative 

The above equations yield 

= [f[{x) + ujf2{x)]gi^^{ux + 6) 

+ [f2i^)+V^fiix)]g2,^iujx + 6). (6) 

The present two-point approximation to ^'{x) which corresponds to the expo- 
nentially fitted formula (2.28) of p| reads 

$2(x) = jj^{[fi{x + h)-Mx-h) + 2Xf2{x)]g,,,{u;x + 6) 

+ [f2{x + h)~ f2{x -h) + 2r]Xf,{x)]g2,,{u;x + 6)} , (7) 

where 

X = huj . (8) 

If the non-oscillatory factors fi{x) and f2{x) are at least three times differentiable, 
then the discretization error associated to (^) shows Olh"^) accuracy 

e'^ix) = -]:hyi^\x + eih)g,,^{ujx + 6) + fi^\x + e2h)g2,rj{uJX + 6)], 
o 

-1< 01,^2 <1. (9) 

Similarly, the present four-point approximation to which corresponds to 

the exponentially fitted formula (3.1) with coefficients (3.3) of reads 

Kix) = ^{[f,{x-2h)-8hix-h)+8h{x + h) 

-fiix + 2h) + 12\f2{x)]gi,^iu;x + 5) 
+ [f2ix - 2h) - 8/2(3; ~h) + 8/2(x + h) 

- f2{x + 2h) + l2T^\h{x)]g2,r,{ujx + 5)} . (10) 

This expression holds for non-oscillatory factors fi{x) and f2{x) which are at least five 
times differentiable and associate an 0{h'^) discretization error 

e'^{x) = ^/i^[/f)(x + ^3%i,,(^x + 5) + /f (x + ^^4%2,,(^X + 5)], 

-1< 03,^^4 <1. (11) 
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2.1.2 Second Order Derivative 

For the second order derivative of (|^), we get 

$"(0;) = [f'!{x) + 2u;f!,{x) + r]u;'Mx)]g,,,{u;x + 6) 

+ [f^{x)+2r]u;f[{x) + r]ij'Ux)]g2,,{u;x + 6). (12) 

hence the approximation to $"(x) which corresponds to the the exponentially fitted 
formula (3.16) with coefficients (3.18) of is provided by the expression 

•^2^ = ^{{[fi{x + h) + {7^X'-2)Mx) + Mx~h)] 

+ A[/2(a; + h) - f2ix- h)]}gi.n{ujx + 5) 
+ {[f2ix + h) + ivX^ - 2)h{x) + h{x - h)] 

+ r/A[/i(x + /i) - /i(x - h)]}g2,r,{oox + 5)]. (13) 

This expression is useful provided the non-oscillatory factors fi{x) and /2(a;) are at 
least four times differentiable. Then its associated discretization error shows 0{h^) 
accuracy 

e'ii^) = -^h^{[f['\x + e,h)+Au:f!i\x)]g^,,{ux + 5) 

+ [f^\x + e,h) + Ar^uff\x)]g2,r,{ujx + 5)} , -l<e,,0e<l. (14) 

The dependence of the discretization errors on the frequency parameter uj will be 
discussed in the case of trigonometric functions {rj = —1) only. Then | cos{ujx + 6)\ < 1 
and I sin(u;x + S)\ < 1, hence Eqs. ([TTl) and ( |T3|) result in the following upper bounds, 
which are independent of uj in the case of Grst order derivatives and hnearly increasing 
with \uj\ in the case of second order derivatives: 

M'^ix) = -h''[\f['\x + d,h)\ + \fi'\x + dM]. (15) 
M',{x) = ^h'[\f['\x + 9,h)\ + \f^'\x + 9,h)\], (16) 

M^{x) = ^hm\fi'\x + e,h)\+4\u;\\f^'\x)\] 

+ [|#(x + WI+4k||/f)(x)|]}. (17) 

2.2 A numerical example 

We consider the case study provided by the test function (2.30) of p[, namely, 

^x) = f{x)cos{ujx) , fix) = 1/(1 + x) . (18) 
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and compute, similar to [Q, the errors associated to the discretized Leibniz and expo- 
nentially fitted approximate derivation formulas at x = 1, under a step-size h = 0.1. 
The computation is done for frequency parameter values a; G [0, 80], using the uniform 
sampling 

ujk = khu;, huj = 0.l. (19) 
In figures 1 and 2 of 0, the linearly scaled errors 

^ r $'(1) - m) < 1 .20) 

I ($'(!)_ otherwise, ci = 2,4, ^ ' 

have been plotted. To allow straightforward comparison of the two methods, the u- 
dependence of the scaled errors (EDI) are shown in figures IT] and 13 below. 



To demonstrate the uniform bounds (^) and (^) respectively, the absolute errors 
associated to the 0{hf') {d = 2, 4) Leibniz formulas (^ and (p!0|), 

E', = $'(1) - m) , (21) 

are also plotted in figures [l| and ^. 

There are three characteristic features which follow from the obtained data. 



• The absolute errors associated to the Leibniz derivation formulas (0) and (|T0|) 
are finite everywhere, irrespective of the value of the frequency parameter u, in 
agreement with the bounds (|15[) and (|T6|). 

The magnitudes of the absolute errors associated to the exponentially fitted 
derivation formulas of reference get arbitrarily large in the neighbourhood 
of the critical points = /cTT , /c = ±l, ±2, ... 

• Let 

uJk = uJo + 2k7i, k = 0, ±1, ±2, ••• (22) 

For each of the two Leibniz approximate derivation formulas, Eqs. (|^ and (|T0|) 
respectively, Eqs. (^) and ( PTj ) yield 

E',M = E'.iuo) , (23) 

i. e., the absolute errors (|2l| ) remain the same over the set of equally spaced 
values (0). These errors show therefore periodic behaviours, of constant specific 
amplitudes, with respect to the frequency parameter uj. Thus, in Fig. |I], the 
amplitude 0.627 x 10~^ equals the accuracy of the two-point derivation formula of 
the factor f{x) in the test function (|T^), while in Fig. |^, the amplitude 0.633x10"^ 
equals the accuracy of the four-point derivation formula of the same factor f{x). 
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Figure 1: The u) dependence of the accuracies of the Leibniz derivation formula (|^) 
and of the two-point exponentially fitted formula (2.28) of 0. The linearly scaled 
absolute errors of (|^) (solid line) damp out as u increases. The linearly scaled errors of 
the exponentially fitted formula (linespoints) are small over narrow ranges of u values 
only. The genuine absolute errors (PT| ) of the Leibniz derivation formula (|^) (dashed 
line) demonstrate the uniform bound (p!5|). 

By contrast, the magnitudes of the absolute errors associated to the exponen- 
tially fitted formulas of the first order derivatives linearly increase with \uj\ over 
the set (^). Hence the accuracy of these exponentially fitted formulas linearly 
deteriorates with the increase of Ic^l. 

• If in Eqs. (^) and ([TT|) the values 9i = 02 = 9^ = 64 = 1 are assumed, then 
leading error estimates are obtained. At the scale of the present figures |l| and |, 
the exact errors ( piD are practically indistinguishable from these estimates (com- 
pare the amplitudes A2{exact) = 0.627 x 10~^ with A2{est.) = 0.625 x 10"^, 
and Aiiexact) = 0.633 x 10"^ with A^iest.) = 0.625 x 10"^). This is a conse- 
quence of the particular test function (|^) proposed in and it should not be 
overemphasized. 
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Figure 2: The uj dependence of the accuracies of the Leibniz derivation formula ( [T0| ) 
and of the four-point exponentially fitted formula (3.1), (3.3) of The linearly scaled 
absolute errors of ( pUD (solid line) damp out as uj increases. Over the sampling (|T^), 
the linearly scaled errors of the exponentially fitted formula (linespoints) fall within 
the specified y-range (and are hence visible) at low uj values only. The genuine absolute 
errors (pT]) of the Leibniz derivation formula (|l^) (dashed line) demonstrate the uniform 
bound iWm. 



Figure 1 of [Q, however, shows that, for the same case study, the estimated errors 
of the exponentially fitted formulas may vary significantly from the exact error 
values. 



The bound ( |T7| ) shows that the appropriate quantities for the study of the depen- 
dence of the accuracy of the discretization errors (O) versus uj are the linearly scaled 



errors 



$"(1) - $2(1 



\uj\ < 1 
— $2(1))/^ otherwise 



(24) 



In figure 3 of [0], however, the dependence versus uj of the accuracy of the ex- 
ponentially fitted three-point second order derivative has been plotted in terms of the 
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quadratically scaled errors 



A" 

^2 




M < 1 
otherwise 



(25) 



Figure |^ summarizes the u dependence of the quadratically scaled errors of the 
exponentially fitted and discrete Leibniz derivatives as well as the linearly scaled errors 
of the discrete Leibniz derivative (|T3|). 
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Figure 3: The u! dependence of the accuracies of the Leibniz derivation formula ([l^) 
and of the three-point exponentially fitted formula (3.16), (3.18) of reference [0]. The 
quadratically scaled absolute errors of (|T^) (linespoints) damp out as u increases. The 
quadratically scaled errors of the exponentially fitted formula (solid line) are shown 
for comparison. The linearly scaled absolute errors (|2^) of the Leibniz derivation 
formula (|13]) (dashed line) show constant amplitude with lu at u > 1. 

From these data, we draw the following conclusions on the approximation of the 
second order derivatives within the two approaches. 

• Similar to the result obtained for the first order derivatives, the errors associated 
to the Leibniz formula are finite everywhere, while those associated to the expo- 
nentially fitted formula diverge at the special point set a; = fcvr , k = ±l, ±2, . . . 
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• At frequency parameter values \uj\ > 1, the Leibniz formula is almost everywhere 
better than the exponentially fitted one. The set of arguments ux + S at which 
the accuracy of the latter formula is comparable to that of the former gets mono- 
tonically smaller as the magnitude of u increases, tending to a countable manifold 
at asymptotically large |ci;|. 

• For the present case study, a phase difference, roughly equal to 7r/2, is noticeable 
between the cj-dependences of the errors associated to the two approaches. This 
feature stems from the prevalence of different asymptotic terms in the discretiza- 
tion errors: a term proportional to u'^ cos(a;) in the exponentially fitted formula 
of 0, and a term proportional to ujsin{uj) in the Leibniz formula. 

The evidence accumulated for the first and second order derivatives raises the 
question on the rate of deterioration of the accuracy of the discretization of the n-th 
order derivative of the function (|^) under asymptotically large lj values. 

Within the exponential fitting, this rate gets proportional to therefore this 
approach yields unsatisfactory results for derivatives of any order n, even n = 1, which 
corresponds to the first order derivatives. 

Let fd^\x) denote the 0{h'^) difference formula of the m-th order derivative of 
the regular factor f{x) entering the reference function ([^). Then the discretization 
error associated to the 0{h'^) approximation of equation @ reads 

Ef\x) = - (^) = E ( ^ ) ^1/^"-'^^) - ft'\^U':^it) , (26) 

therefore it behaves like Icu]""^ at asymptotically large lv. 

3 Quadratures 

3.1 Interpolatory Simpson quadrature sums 
for oscillatory functions 

Let 

i::'hb]f= f' f{x)gs,,{uJX + S)dr, s = l,2, r/ = ±l, (27) 

J a 

denote the distinct integrals corresponding to each of the four basic integrands (|l|). 

In this section, we discuss interpolatory Simpson quadrature sums for oscillatory 
functions. These provide an alternative to the exponentially fitted formulas derived 
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in section 4 of 0. We start with the derivation of rehable quadrature sums for the 
solution of the basic integrals ( pTl) . Simpson quadrature sums for integrals 

rb 

/[a, 6]$ = / $(x)cfc, (28) 



with integrands of the form will then be given by linear combinations of the 
basic quadrature sums. 

An important preliminary operation is the mapping of the integration range [a, b] 
onto the fundamental range [—1,1] via the transform of independent variable, 

x = c + hy; ye [-1,1]] c= {b + a)/2; h = {b-a)/2, (29) 

where c and h denote the centre and respectively the half-length of the integration 
range [a,b]. 

As a result of this operation, we get the following relationships which are best 
written in matrix form 

Thus, the calculation of each of the two integrals (^) occurring in the left hand side 
has been reduced to the evaluation of the pair of fundamental integrals 

I^,,[a,b]iP = j\iy)gs,,iXy)dy, ^Piy) = /{c + hy) ; s = l,2, (31) 

with a transfer matrix R of elements 

Rs,^ = hgs,n{'fc) , s = l,2;?7 = ±l. (32) 

Here, A is the dimensionless parameter defined by Eq. (||), while ipc = ujc + 6 denotes 
the phase of the center c of the integration range [a, b]. 

To derive Simpson quadrature sums for each of the basic integrals (0), the 
regular part f{x) of the product (|ID is interpolated by an arc of parabola, L2{x), at 
the knots 

a = Xq < X2 < Xi = b . (33) 

We use the Newton representation of the unique Lagrange interpolation polynomial of 
f{x) at the knots (||) ||. 

The generalization of the classical Simpson quadrature to the solution of the 
integral (^) will then be 

Q^;;[a,b]f = i::^'[a,b]L2= f L2{x)gs,^{ujx + 5)dx . (34) 

J a 
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Taking into account (|30|), the computation of the last integral is reduced to that 
of the pair of integrals J^^[— 1, l]i2 and l2^^[—l, 1]£2, where ^2(1/) = -^2(c+ hy). To get 
well-conditioned expressions of these integrals, £2(1/) is expressed in terms of Chebyshev 
polynomials of the first kind, 

^2(l/) = EA2T,(i/), (35) 

1=0 

with the coefficients given respectively by 

(3o2 = |[(3 - l/p)/o + (2 + p + l/p)/2 + (3 - p)/i] , 
P12 = K/1-/0), 

(522 = + l/p)/o - (2 + p + l/p)/2 + (1 + p)/i] . (36) 
fk = f{xk), fc = 0,l,2 (37) 



Here 
while 



1 + ^/2 

, 38 

1-1/2 

denotes the ratio of the lengths of the two subranges created by the inner abscissa 
1/2 = (x2 — c)/h inside the fundamental range [-1, 1]. 

General expressions of the reduced integrals J^^[— 1, 1]£„ have been derived in 
[There, however, two misprints have to be corrected. First, the right hand side of the 
expression (28) of the coefficients Sq^2k-i is to be preceded by a minus sign. Second, 
in Eq. (39), the summation is to proceed from 1 to k, over hypergeometric functions 
o-Fi(j + I; ^rjX"^).] In view of the low polynomial degree of the expansion (pS]), equa- 
tions (38) and (39) of 0] provide suitable expressions for the moments of the involved 
Chebyshev polynomials. 

The obtained expression of J^^[— 1, l]i2 is transformed to a well-conditioned one 
by use of the recurrence relation (40) of and we finally get 

Ji\[-1, l]i2 = C,{p) ■ oFid; ir/A^) + 4^^'(/o + /i)oi^i(i; i^A^) , (39) 
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IY2 = lx{fi - /o)oFi(f ; ir/A^) , (40) 



where 



dp) = l\{2- l/p)/o + (2 + p + l/p)/2 + (2 - p)/il . (41) 
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If the knots (p3D are equally spaced, then p = 1 and the coefficient (^l]) goes into the 
classical Simpson quadrature sum for the regular factor /(x), 



Ci(l) = -(/o + 4/i + /2). (42) 

Similar to the interpolatory quadrature sums at Clenshaw-Curtis and related 
abscissas, [Q, the only practical difficulty in the implementation of these quadrature 
sums into a code is the accurate computation of the hypergeometric functions o-^i of 
interest. The preparation of a paper which describes the computer programme devoted 
to the computation of basis sets of hypergeometric functions o-^i to machine accuracy 
is underway and it will be submitted to this journal in the nearest future. 

Equations ([3^), (|30D, (^), (|40|), (^) provide the required generalization of the 
classical Simpson quadrature rule to each of the four reference functions (|l|). Indeed, in 
the limit u; — > 0, these equations result in Q^]fj[a, b]f = |/i(/o + 4/2 + fi)gs,r){S). On the 
other side, in the same limit, Eq. (|l]) yields F(x) = f{x)gs,-q{S), hence Fk = fkgs,ri{S)- 
Therefore, Q^i^[a, b]f = \h{FQ + 4^2 + Fi), which is nothing but the classical Simpson 
quadrature sum Q[a,b]F. 

The Simpson quadrature sum for the integral ( p8D with the integrand $(x) given 
by (^ is then given by 

[a,6]<l> = (/^,[-l,l]f2 + /2',hl, 1^2)^1,. 

+ lY'i + vlU-l, 1%)R2,, , (43) 

where d^iy) and ^2(z/) denote the second degree interpolatory polynomials associated 
to the regular factors /i(c + hy) and /2(c + hy) respectively. 



3.2 Leading error estimate 



The error associated to the quadrature sum (Q), which replaces the computation of 
the basic integral (^Tf ), is given by 

/,^;f[a,6]A/ = j' Af{x)gs,,{u;x + 6)dx, 

A/(x) = f\x)-L2{x). (44) 



In view of the relationship (|30|), an estimate of this expression requires two esti- 
mates (pll), s = 1, 2, for the function 

A^(y) = ij{y) - i,iy) = fic + hy) - ^v)- (45) 
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To get the leading terms of the two estimates, the reference function /(c + hy) 
is expanded in Taylor series up to fourth order around the point x = c while the 
coefficients of the interpolation polynomial are expressed in terms of the 0{h'^) ap- 
proximate derivatives /2(c) and /2(c) of f{x). After straightforward algebra, we get 
the Chebyshev polynomial expansion 

/\^{y) ^ \a,[T2{y) - T^{y)] + \a^[T^{y) - T^{y)] + \a:,[TM - T,{y)] , (46) 
where 

«o = ^[(1 - l/p)/o + (2 + p + l/p)/2 + (1 - p)/i] - /(c) 

~ -^(3 + p)y2a3-l/2«4, (47) 

«3 = \h'f'\c), (48) 

«4 = ^/^V^'^c). (49) 

The coefficient vanishes identically in the case of equally spaced knots. It arises 
only within the generahzed Simpson quadrature sum (|3^) under p 7^ 1. 

The leading error estimates associated to the quadrature sums (|39D and (|40|) are 
now immediate: 

/i\[-l,l]A^ = -^[(ao + «4)-oi^i(f;ir7A=^)-^«4-oi^i(i;ir/A2)], (50) 
/2\[-l, = -^\a, ■ oFi(i; ir^A^) . (51) 

Taking into account the expressions (l47|) , (^Sf ), and (|49|) , we conclude that the 
equally spaced Simpson quadrature sums ( p9D and (^) show 0{h'^) accuracy. However, 
under a nonuniform mesh, the order of accuracy is reduced to 0{h^), due to the non- 
vanishing ao coefficient. 

Let in (|50|) and (|5lD the hypergeometric functions o-^i(f ; jV^"^) o-^i(|; jV^"^) 



be expressed in terms of the fundamental hypergeometric functions 0-^1 (|; jV^^) (which 
equals cos(A) under r] = —1, respectively cosh(A) under t] = 1) and o-Fi(|; i^^A^) (which 
equals sin(A)/A under t] = —1, respectively sinh(A)/A under rj = 1). These operations 
put into evidence the occurrence of an u!~^ power law at asymptotically large u, com- 
ing from the hypergeometric functions. Therefore, the accuracy of the interpolatory 
Simpson quadrature sums is expected to increase with increasing u. 
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3.3 An accuracy test 



A meaningful comparison of the above interpolatory Simpson and exponentially fitted 
quadrature sums of is provided by the test case (4.18) of reference 0, that is, the 
integral (El) of the oscillatory function 



$(a:) = — /^(x) cos(co'x) — ujf{x) sin(u;x), /(x) = 1/(1 + x) , 

over the range [a, b] = [0.9, 1.1]. The primitive of the integrand ( p^ is the function 
hence the exact value of I[a, 6]$, Eq. (pSf), is immediate: 



(52) 



J[a,6]<l> 



rh 1 + c 1 
- cos(v9c) cos(A) H — sm{ipc) sin(A)J / 



1 + c 



(53) 



where c = 1 ; h = 0.1 . With the quadrature sum Q[a, 6]$ given by (^31), we can compute 
the absolute quadrature errors, at various values of the frequency parameter u, 



A = /[a,6]$-g[a,6]$, 



(54) 



An important feature stemming from equation (^) is the occurrence of two 
distinct periodic oscillatory behaviours of the considered integral. First, there is a 
characteristic (pc-dependent periodicity, defined by the phase (pc = ujc + 6 of the centre 
c of the integration range through the factors cos{(pc) and sm{(pc)- Here, this kind 
of periodicity is characterized by short wavelength oscillations of period T^^ = 27r. 
Second, there is a characteristic X-dependent periodicity, defined by the phase A of the 
integration range half-length through the factors cos(A) and sin(A). Here, this kind of 
periodicity is characterized by long wavelength oscillations, of period Tx = 2tt / h = 207r, 
which is ten times larger than T^p^. 

Figure |^ plots the obtained dependence A versus oj over the frequency parameter 
range u G [0, 500]. To allow straightforward comparison with the optimal exponentially 
fitted P = 1 method of 0, a supplementary plot concerning the uj dependence of the 
quadrature errors of this method is given in figure |^. In both figures, the sampling (|T9|) 
was used. Separate plots have been needed simply because the relevant y-scale of 
figure ^, (—8.0 x 10~^, +8.0 x 10~^), is two orders of magnitude smaller than that 
required in figure ||, (—1.0 x 10~^, +1.0 x 10"^). 

The inspection of the obtained numerical evidence results in the following conclu- 
sions concerning the interpolatory Simpson and exponentially fitted quadrature sums. 

• Both kinds of periodicities discussed above are present in figures ^ and ^. While 
the yjc-induced periodicity pattern is the same in the two figures, the A-induced 
periodicity pattern is substantially different. 
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Figure 4: Absolute quadrature errors, A, Eq. (|5lD, associated to the interpolatory 
Simpson method for the oscillatory function (^) (solid line). The two envelopatrices 
A = ±0.0025/u; (dashed lines) illustrate the asymptotic error behaviour established by 
the error analysis. 

• Being controlled by the function rjoi^—X^) = sin(A)/A, the A-periodicity of the 
exponentially fitted method shows, over sets of points uj > Tx/2, a period exactly 
equal to Tx- 

The error magnitudes are roughly the same over sets of u values separated by 
multiples of Tx. 

• By contrast, the A-periodicity of the interpolatory Simpson quadrature sums is 
controlled by higher order hypergeometric functions o-^i, the periods of which 
only asymptotically tend towards the value Tx- (Thus, for the present problem, 
approximately periodic long wavelength oscillations can be defined starting with 
values u > Tx- As a; increases, the quasi-period lengths monotonically decrease 
towards Tx- At values u > 400 an agreement with Tx within values better than 
one percent is obtained.) 

The behaviour of the error amplitudes with ui is in agreement with the error 
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Figure 5: Absolute quadrature errors associated to the optimal exponentially fitted 
P = 1 method of 0. 



analysis done in subsection |3.2| . The important detail to be noted here is the 
occurrence, in the integrand of a non-oscillatory term which is proportional 
to uj. This results into an asymptotic uj~^ power law behaviour of the A-defined 
long wavelength amplitudes. In figure the occurrence of this asymptotic regime 
is visualized by means of the two envelopatrices A = ±0.0025/ci;, which provide 
true upper bounds to the discretization errors at u values beyond T^. 

In reference , breakdown of the quadrature errors was noticed to occur around 
A = 2k7T for the suboptimally fitted method P = and around the critical values 
X = kn {k = 1, 2, ...), for the optimally fitted method P = 1 (only the latter 
being visible at the error scale, —0.01 < Ap < 0.01, of fig. 4 of reference 0). 
Figure ^ is consistent with this frame. An increase of the mesh density ( [TP| ) adds 
supplementary details in the neighbourhood of the critical points only. 
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4 Conclusions 



The approximate p-point Leibniz derivation formulas and the interpolatory Simpson 
quadrature sums discussed in this paper provide simple, versatile, and efficient tools 
for the analysis of oscillatory phenomena. 

Both theoretical considerations and numerical evidence concerning the depen- 
dence of the discretization errors on the frequency parameter of the oscillatory func- 
tions show that the accuracy gain of the present formulas over those based on the 
exponential fitting is overwhelming: 

1. The discretization errors of the present numerical differentiation and integration 
formulas are finite everywhere with respect to the frequency parameter u, whereas 
those of the corresponding exponentially fitted ones diverge at specific countable 
sets of uj values. 

2. Even if the huge errors coming from convenient neighbourhoods around diverging 
points of the exponentially fitted methods are disregarded, the accuracies of the 
present methods still remain, in the average, two orders of magnitude better with 
respect to those based on the exponential fitting. 

3. Over sets of u values separated by entire periods of the oscillatory function, 
the absolute errors of the approximate Leibniz n-th order derivatives increase 
following an |ci;|"~^ power law, as compared to the la;!"" power law specific to the 
exponentially fitted formulas. 

Thus, the Leibniz first order derivatives (|^) and (|TUp are characterized by the 
uniform bounds (^) and ( |T6| ) respectively, whereas the exponentially fitted ones 
linearly deteriorate with uj. As it concerns the second order derivative, the Leibniz 
formula ( [T3| ) shows linear deterioration with u, Eq. (|1^, as compared to the 
quadratic deterioration with uj of the exponentially fitted counterpart. 

4. The errors associated to the interpolatory Simpson quadrature sums ( ^Hj ) show 
quasi-periodic behaviour with respect to uj, with a damping uj~^ power law of the 
amplitudes, coming from the hypergeometric functions o-F^i- As a consequence, 
the accuracy improves as u increases over sets of values separated by hypergeo- 
metric function induced quasi-periods. 

The best exponentially fitted Simpson quadrature formula (P= 1) of reference 
yields comparatively large and roughly uniform error magnitudes over sets of uj 
values separated by multiples of the A-induced period Tx. 
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Of course, the exponential fitting provides approximating formulas are signifi- 
cantly better in comparison with classical formulas (wherein the oscillatory character 
of the function of interest is ignored). However, in view of the abovementioned results, 
the claims made in section 6 of reference that the exponentially fitted formulas are 
"working optimally" on functions of the form (j^) and that further useful extensions 
can be developed concerning the numerical differentiation and integration, are to be 
regarded with caution. 

We end with an interesting observation concerning the interpolatory Simpson 
quadrature sums. They can be optimally extended to five-point interpolatory quadra- 
ture sums to yield quadrature rules along the lines of QUADPACK |^ . These can be then 
conveniently implemented in automatic quadrature codes able to match the advan- 
tages of the polynomial function approximation at Clenshaw-Curtis or Gauss-Kronrod 
quadrature knots with the use of the whole information acquired on the integrand at 
previous stages of the adaptive subrange subdivision. This point will be discussed 
separately in a future paper. 
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